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We discuss the implementation of a directed geometrical worm algorithm for the study of quantum 
link-current models. In this algorithm Monte Carlo updates are made through the biased reptation 
of a worm through the lattice. A directed algorithm is an algorithm where, during the construction 
of the worm, the probability for erasing the immediately preceding part of the worm, when adding 
a new part, is minimal. We introduce a simple numerical procedure for minimizing this probability. 
The procedure only depends on appropriately defined local probabilities and should be generally ap- 
plicable. Furthermore we show how correlation functions, C(r, r) can be straightforwardly obtained 
from the probability of a worm to reach a site (r, r) away from its starting point independent of 
whether or not a directed version of the algorithm is used. Detailed analytical proofs of the validity 
of the Monte Carlo algorithms are presented for both the directed and un-directed geometrical worm 
algorithms. Results for auto-correlation times and Green functions are presented for the quantum 
rotor model. 

PACS numbers: 02.70.Ss, O2.70.Tt, 74.20.Mn 
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I. INTRODUCTION 



Improving and developing new numerical algorithms 
lies at the heart of computational physics. Amongst 
others, Monte Carlo (MC) methods are often seen as 
the best choice for the study of phase transitions tak- 
ing place in classical or quantum models. For the study 
of spin models for example, cluster algorithms, either in 
the classical 0,0 or quantum [J HJ case, perform non- 
local moves in phase space, allowing for the treatment of 
systems much larger than with traditional local update 
methods (single spin- flip algorithms). These types of al- 
gorithms have almost completely solved the problem of 
critical slowing down arising near phase transitions. 

The class of systems for which cluster methods are 
known to exist is limited, and it is therefore of great inter- 
est to search for new algorithms possessing the same effi- 
cient features for other models. In this context, we have 
proposed recently a non-local "worm" algorithm for the 
study of quantum link-current models |fj . These models 
arise from a phase approximation of bosonic Hubbard 
models, but are also relevant in the context of quan- 
tum electrodynamics 0. Previous MC simulations of 
the quantum link-current (quantum rotor) model used a 
local algorithm suffering from critical slowing down. In 
the new algorithm |(j updates are made by reptating a 
"worm" through the lattice 0,0- Since the movement 
of the worm only depends on a few probabilities deter- 
mined locally with respect to the current position of the 
"head" of the worm, we call this type of algorithm a ge- 
ometrical worm algorithm as opposed to other recently 
developed worm algorithms based on high-temperature 
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series expansions [||. The geometrical worm algorithm 
gives rise to very small autocorrelation times and by di- 
recting the algorithm these autocorrelation times can be 
even further reduced. 

In this paper, we briefly recall the principles of the 
geometrical worm algorithm During the construc- 
tion of a worm a new part is added to the worm by 
moving the worm through one of the a nearest neigh- 
bor links. Usually the associated a probabilities, p a , 
are chosen in an un-biased geometrical way and there 
is therefore a significant probability that the new part of 
the worm will back-track in its own path, thereby eras- 
ing the immediately preceding part. In many cases this 
back-tracking (or bounce) probability is the dominant 
probability among the a probabilities and using these 
un-biased probabilities is therefore clearly rather waste- 
ful. Here we describe an improvement of this geometric 
worm algorithm, which we call the directed worm algo- 
rithm, as a reference to recently developed directed loop 
methods for Quantum Monte Carlo simulations of spin 
systems |9j- This directed geometrical worm algorithm 
is identical to its un-directed counterpart except for the 
fact that the probabilities p G are now chosen in a biased 
way, using knowledge of the immediately preceding step 
in the construction of the worm. These biased proba- 
bilities can all be tabulated at the start of the simula- 
tion and the additional computational effort stems solely 
from the significantly wider distribution of the directed 
worms. The directed algorithm gives rise to even better 
results, as will be shown in the following part of this pa- 
per, where we present results on autocorrelation times for 
both directed and "undirected" worm algorithms. The 
procedure for choosing the "biased" probabilities lead- 
ing to the directed algorithm is quite general and should 
be applicable to other algorithms that depend on local 
probabilities. Furthermore, we show how Green func- 
tions C(r, r) of the original quantum model can be mea- 
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sured efficiently during the construction of the worm by 
calculating the probability that the worm reaches a given 
site (r, r) away from its starting point, independently of 
whether a directed or un-directed algorithm is used. For 
both the derivation of the directed algorithm and the 
measurements of correlation functions, analytical proofs 
of the validity of the algorithms are presented. 

The outline of the paper is as follows: in the next sec- 
tion, we present the quantum rotor model and introduce 
some useful notation. Then, a brief description of the 
" undirected" worm algorithm is given in section IIII Al 
before proceeding to the main contents of this paper, a 
description of the directed geometrical worm algorithm 
A simple procedure for numerically de- 
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termining the biased probabilities for the worm moves 
is presented. In addition, we derive a proof of detailed 
balance for the directed worm algorithm. In order to 
compare our algorithms to related ones, we present in 
section IIII CI another recent approach due to Prokof 'ev 
and Svistunov originally based on high temperature 
series expansion for classical statistical models, which we 
therefore will refer to as "classical worms" throughout 
this paper. In section llVl we estimate the efficiency of 
the algorithms by calculating auto-correlation times in 
the MC simulation, and compare to both undirected and 
classical algorithms. We then discuss the measurements 
of correlation functions within the worm algorithm in sec- 
tion and show some results at a specific point of the 
phase diagram. We conclude with a discussion of the 
features of the directed algorithm. 



II. THE MODEL 

Many magnetic systems, Josephson Junction arrays 
and several other systems can be described by a quantum 
rotor model yjfl : 



1 d V d x 



cos W r — t 
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Here, 9 r is the phase of the quantum rotor, t the renor- 
malized coupling strength and \x v an effective chemical 
potential. If /i = it can be shown that this model 
displays the same critical behavior as the D + 1 dimen- 
sional XY-mode\. However, when fi r ^ this model is 
not amenable to direct numerical treatment in this rep- 
resentation due to the resulting imaginary term. It is 
therefore very noteworthy that an equivalent completely 
real representation in terms of link-currents exists even 
for non-zero fi r . This link-current (Villain) representa- 
tion is a classical (2+l)D equivalent Hamiltonian that is 
usually written in the following manner [Tlj : 
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The sum is taken over all divergenceless current configu- 
rations V • J = 0. The degrees of freedom are "currents" 
J = (J x , J v , J T ) living on the links of the lattice. These 
link variables J x , J v , J T = 0, ±1, ±2, ±3 . . . are integers. 
K is the effective temperature, varying like t/U in the 
quantum rotor model. We refer to Ref. 11] for a pre- 
cise derivation of this model and for a description of its 
physical implications. 

Another incentive for studying the critical behavior 
of the quantum rotor model comes from the close re- 
lation between this model and bosonic systems. Bosonic 
systems with strong correlations are often described in 
terms of the (disordered) boson Hubbard model: 7?bH = 
E r (f "r - A*r" r ) - t E(r,r') (H®* + c.c) . The corre- 
lations are here described by U the on-site repulsion. 
The hopping strength is given by to and fi r the chem- 
ical potential varying uniformly in space between /i ± A. 
h T = &l$r is the number operator. If we set <3E> r = 
|$ r |e ier and integrate out amplitude fluctuations, it can 
be shown that -ffbH is equivalent to the quantum rotor 
model For systems where amplitude fluctuations 

can be neglected at the critical point, such as granu- 
lar superconductors and Josephson junctions arrays, the 
quantum rotor model should therefore correctly describe 
the underlying quantum critical phenomena. 

In the following we only discuss the quantum rotor 
model in d = 2 dimensions corresponding to the d + 1 
dimensional link-current model. In general a non-zero fi 
will introduce separate dynamics for the time and space 
directions from which a dynamical critical exponent, z 
can be defined. If the divergence of the spatial correlation 
length close to criticality is characterized by the exponent 
v, z is defined by requiring that the correlation length in 
the time direction diverges with the exponent zv. For 
what we will be discussing here (i = and z = 1. 



III. ALGORITHMS 

A. The Geometrical (Undirected) Worm algorithm 

The quantum rotor model has been extensively stud- 
ied in the link-current representation using conventional 
Monte Carlo technique using local updates 0, 0, 0, 

mm 

Conventional Monte Carlo updates on the model (J2J) 
consists of updating simultaneously four link variables 
as shown in the left part of Fig. ^ (A). To ensure er- 
godicity, one also has to use global moves, updating a 
whole line of link variables (B in Fig. The accep- 
tance ratio for these global moves becomes exponentially 
small with the system size for large systems. Many inter- 
esting quantities such as the stiffness, necessary for the 
determination of the critical point, and current-current 
correlations, necessary for the calculation of transport 
properties such as the resistivity and the compressibility, 
are only non-zero when these global moves are successful. 
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FIG. 1: Monte Carlo moves in the model. To ensure the 
divergenceless condition, only closed moves can be performed. 
On the left part of the figures, previous Monte Carlo updates 
with the local algorithm are depicted. On the right part is 
given an example of a move with the worm algorithm starting 
from an initial random site (black dot). 



An effective Monte Carlo sampling of these global moves 
is therefore imperative and it is easy to understand that 
the performance of the local algorithm is rather poor, es- 
pecially near a phase transition: critical slowing down in 
the Monte Carlo simulations prohibits the study of large 
system sizes. 

In order to be able to correctly describe the directed ge- 
ometrical worm algorithm we have to review the geomet- 
rical worm-cluster algorithm introduced in reference [|J 
in some detail. This algorithm allows for non-local moves 
as the ones depicted on the right part of Fig. (C). The 
performances of this algorithm have been reported in the 
previous work [6|. This algorithm is closely related to 
other cluster algorithms 0,0)01 an d especially to "worm" 
algorithms 0, 0, @, from which we have borrowed the 
name. We stress that this algorithm is different from 
the "classical" worm algorithm presented in Ref. [8( in 
the sense that it is geometrical: link variables are not 
"flipped" with a thermodynamic probability, instead, a 
new part of the worm is added by selecting a direction 
according to a locally determined probabilities p a . Since 
these probabilities only depend on the local environment 
we call them "geometrical' probabilities. Secondly, even 
though the local probabilities p„ do depend on the ef- 
fective temperature, K, we always have J2aP<? = 1 (P er 
definition) and the only effect of the temperature is there- 
fore to preferentially move the worm in one direction as 
opposed to another one. In some sense this is very sim- 
ilar to the N-fold way of performing Monte Carlo 
simulations. 

We now describe the contents of the geometrical algo- 
rithm: we update the configurations by moving a "worm" 
through the lattice of links. The links through which 
the worm pass are updated during its construction. The 
configurations generated during the construction ("rep- 
tation" ) of the worm are not valid (the divergenceless of 
J is not fulfilled) but at the end of its path, when the 
worm forms a closed loop, this condition is verified and 



the final configuration is valid. 

We first define a convention for the orientation of the 
lattice. Around each site with coordinates (r = (x, y),r), 
there are six links on which the integer currents JZ T s are 
defined with a — x,y, r, —x, —y, —r. The change in JZ ^ 
that the worm will perform during its course depends on 
whether a is an incoming or outgoing link: here our con- 
vention is to consider positive x,y,T as outgoing direc- 
tions and —x, —y, —t as incoming (see left lower part of 
Fig. P. 

If the worm is leaving the site (r, r) passing through 
an outgoing link a = X,y, r, then 

If it is leaving through an incoming link a = —x, —y, — r 
, we have 

J(r,r) "> J(r,r) ^ (4) 

Graphically, the convention means that the update in 
JZ, r \ is +1 (—1) if the worm goes in the same (opposite) 
direction as the arrows denoted in the left lower part of 
figure 2] 

The construction ("reptation") of the worm can be de- 
scribed the following way: first, we start the worm at 
a random site Ni = (ri,ri) of the lattice (black dot in 
figure [Q. From this site, the worm has the possibility to 
go to one of its six neighboring sites. To choose which 
direction to take, a weight A a is calculated for the six 
directions a = ±x, ±y, ±r. For A" , we use a Metropolis- 
like weight: 

A' =min(l 1 exp(-(££ (5) 

where E a = ~ (J CT ) 2 — [iJ° '8 a , T is the local energy carried 
by the link a and E' a = ±( J a ± l) 2 - fj,( J a ± 1)£ CT>T is the 
local energy on the link a if the worm passes through this 
link. The plus or minus sign depends on the incoming 
or outcoming nature of the link (see above). Please note 
that there are other possible choices for A CT [l^ ■ 

Once the A CT 's are calculated, one computes the prob- 
abilities p a by normalizing the weights A,, : 



where N — A a is the normalization. A random num- 
ber uniformly distributed in [0, 1] is generated and a di- 
rection a chosen according to equation J^J. Once a di- 
rection is chosen, the corresponding link variable J CT is 
updated by ±1 and the worm moved to the next lattice 
site in this direction. 

From there, we apply the same procedure to choose 
another site, modify the link variable, move the worm 
until the worm eventually reaches its starting point and 
forms a closed loop. This is then the end of this non-local 
move. 

To satisfy the detailed balance condition, this worm 
move must either be accepted or rejected. To check this, 
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one has to store the initial and final normalizations iV] 
and iV/j (calculated as in Eq. 0) of the weights at the site 
s\ = (ri,Ti). Nl is the initial normalization before the 
worm is inserted and JV/ the final normalization after 
the worm reaches the initial point. The worm move is 
then accepted with probability -W* /iV/ . If the move is 
rejected, we have to cancel all changes of the link-currents 
made during the construction of the worm. During a 
typical simulation the rejection probability is usually very 
small. 

As already mentioned, the link configurations gener- 
ated during the worm move do not satisfy the divergence- 
less constraint, but it is easy to see that the final config- 
uration does. It is important to note that the worm may 
pass many times though the same link and that at each 
step, it can bounce back (back-track) to the previous lat- 
tice site in its path. 

A proof of detailed balance for this algorithm is ob- 
tained by considering the moves of the worm and of an 
anti-worm, going exactly in the opposite direction 
This worm algorithm satisfies ergodicity since the worm 
can make local loops and line moves as in the local algo- 
rithm, which is ergodic. 

All in all, the geometrical un-directed worm algorithm 
can be summarized using the following pseudo-algorithm: 

1. Choose a random initial site si = (ri,Ti) in the 
space-time lattice. 

2. For each of the directions a = ±x, ±y, ±r, 
calculate the weights A a s with A a s = 
mm(l,exp(-AE°jK)), AE%. '= E'f. - E°.' 

3. Calculate the normalization N Si = A°. and the 
associated probabilities p a s . = A°. / N Si . 

4. According to the probabilities, p" s . , choose a direc- 
tion a. 

5. Update the Jf. for the direction chosen and move 
the worm to the new lattice site Sj+i. 

6. If Si ^ s\ goto|21 

7. Calculate the normalizations N Sl and N Sl of the 
initial site, si, with and without the worm present. 
Erase the worm with probability P e = 1 — 
min(l,iV Sl /7V Sl ). 



B. The Directed Geometrical Worm algorithm 

The above algorithm for geometrical worms is not op- 
timal since the worm quite often will choose to erase 
itself by returning to the previous site. While it is in 
general not possible to always set this back-tracking (or 
bounce) probability to zero it is quite straightforward to 
choose the probabilities p a s . such that the bounce or back- 
tracking probability will be eliminated in almost all cases 



and in general will be as small as possible. The proce- 
dure for doing this amounts to solving a simple linear 
programming optimizing problem. If we consider models 
with disorder this has to be done at each site, but the 
correctly optimized (biased) probabilities p". can still be 
tabulated at the start of the calculation. 

In order to see how we can minimize the back-tracking 
probability let us define the 6x6 matrix P Si of proba- 
bilities where the element P s . of the matrix P Si is given 
by the conditional probability p Si {crk\vi) for going in the 
direction o~k at the site Sj if the worm is coming from the 
direction o~i. The back-tracking probabilities at the site 
Si now correspond to the diagonal elements of the matrix 
P Si . For the algorithm described in the previous section 
Psi(o~k\o~i) was simply chosen as A°* /N Si independent of 
a i . Thus all the columns of P Si were the same and P Si 
had in general rather large diagonal elements. However, 
as we shall see below, the matrix P Si only needs to satisfy 
the following two conditions in order to define a working 
geometrical worm algorithm. These conditions are: 

y^.Psj (<Tk\<ri) = 1 (probability) (7) 

k 

p s k ! PsMm) t 

—rr = — — — - = — detailed balance . 8 

"if PsMiWk) a s \ 

These conditions are not very restrictive and will in most 
cases allow us to define a matrix P Si with all the diagonal 
elements (back-tracking probabilities) equal to zero. The 
conventional geometrical worm algorithm, discussed in 
the previous section, corresponds to P^ = A° h /N Si . 

If we define a function / as the sum of the diago- 
nal elements of P Si , f — J2kPs k > we can reformulate 
the search for a matrix P Si with minimal diagonal ele- 
ments as a standard linear programming problem. Writ- 
ing P^ k = 1 — X^/c Pal we should minimize / subject 
to the constraints X)z^/c Pal — 1 The minimum can 
be found using standard techniques of linear program- 
ming and corresponds in almost all cases to / = 0. 
The matrix P Si depends on the value of all the 6 link- 
currents J£. . During the construction of the worm only 
sites where \j~ x + J~* + J-J - J* - Jf. - JJ. = 1 , s, =f si 
occur. Since in general \Jf. \ will almost never exceed a 
certain value J ma x it is easy to construct a lookup table 
for the matrices P Si at the beginning of the simulation 
and only calculate P Si ({ Jf. }) during the simulation if for 
some o~ |JfJ > J max . 

This idea of minimizing the bounce processes is also at 
the heart of Quantum Monte Carlo directed loop meth- 
ods 0. The previous restrictions on the matrix P and 
the way to solve them numerically are indeed very gen- 
eral, and constitute a simple framework for how one can 
construct a directed algorithm out of a "standard" non- 
local loop, worm or cluster algorithm. 

We can now define a directed geometrical worm algo- 
rithm with minimal back-tracking probability. Using a 
pseudo-code notation we have: 
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1. Choose a random initial site si = (ri,n) in the 
space-time lattice. 

2. If i = 1 then: For each of the directions cr = 
±a;, ±y, ±t, calculate the weights A a s . with A". = 
min(l,exp(-AE^./K)), AE'J. = E'jf - E°.. Cal- 
culate the normalization N Si = ^2 a A a . and the 
associated probabilities p° — A° /N Si . Else: Ac- 
cording to the incoming direction, ai , set p" Si equal 
to the Z'th column of P Si . 

3. According to the probabilities, p°., choose a direc- 
tion a . 

4. Update J a . for the direction chosen and move the 
worm to the new lattice site Si+i. 

5. If Si ^ si gotoEl 

6. Calculate the normalizations iV sl , of the site s\ 
with the worm present, and N S1 , without the 
worm. Erase the worm with probability P e — 
l-mm(l,N Sl /N Sl ). 

Now we turn to the proof of detailed balance for the 
directed algorithm. Let us consider the case where the 
worm, w, visits the sites {s± . . . sn} where s\ is the initial 
site. The worm then goes through the corresponding link 
variables {1% . . . In}, with li connecting s, and Sj+i. Note 
that sn is the last site visited before the worm reaches s\. 
Hence, sn and si are connected by the link In- The total 
probability for constructing the worm w is then given by: 




The index a denotes the direction needed to go from s% 
to s 2 , P S1 is the probability for choosing site s 1 as the 
starting point and P^ is the probability for erasing the 
worm w after construction. p Sj (sj_|_i|sj_i) is the condi- 
tional probability for continuing to site S;+i, at site Sj, 
given that the worm is coming from Sj_i. If the worm w 
has been accepted we have to consider the probability for 
reversing the move. That is, we consider the probability 
for constructing an anti-worm w annihilating the worm 
w. We have: 

P a , = P 9l (l - P?)S- f[ p Sz ( Sl -i\s l+ i)- (10) 

iVsl r=N 

Here, the index a denotes the direction needed to go from 
s\ to sn, Note that, in this case the sites are visited in 
the opposite order, si, Sn, ■ ■ ■ , S2. From Eq. (JSJ we have 
that p Si (s l+1 \s l ^ 1 )/p Sz (s i ^ 1 \s l+1 ) = A a *jA? a \. Since, 

KJKi = exp(-A££/J0> i= 1 ■ ■ -N, (11) 
and since P S1 = P Sl , we find: 

P ^ = \ = 3w^- AE ™ ,K) - (12) 



where Ai?Tot is the total energy difference between a con- 
figuration with and without the worm w present. With 
our definition of P e we see that (1— P e (w))/ (I— P e (w)) = 
N Sl /N Sl and it follows that: 

^- = cxp(-AE Tot /K). (13) 

For a worm of length N there are iV starting sites that 
will yield the same final configuration. The above proof 
shows that for each of the starting sites there exists 
an anti-worm, w, such that P w — exjp(—AET t/K)P lB . 
Hence, if we by [i denote the configuration without the 
worm and v the configuration with the worm and fur- 
thermore let P w (si) denote the probability of building 
the worm w starting from site Si, we see that: 

= exp(-ASTot/^). (14) 

Ergodicity is proved the same way as for the undirected 
algorithm as the worm can perform local loops and wind 
around the lattice in any direction, as in the conventional 
algorithm. 



C. The Classical Worm algorithm 

Prokof 'ev and Svistunov have proposed a very ele- 
gant way of performing Monte Carlo simulations on the 
high temperature expansion of classical statistical me- 
chanical models using worm algorithms. In order to dis- 
tinguish between the algorithms we call this algorithm 
the classical worm algorithm. In a recent study [ljj these 
authors have performed simulations on the quantum ro- 
tor model in the link-current representation, Eq. @. Due 
to the divergenceless constraint, the classical worm algo- 
rithm is in this case quite close to the geometrical worm 
algorithm proposed previously in Ref. 6 and not directly 
related to the high temperature expansion of this model. 
Recasting their algorithm in the same framework used 
above we outline our understanding of their algorithm 
below for comparison: 

1. Choose a random initial site si = (ri,ri) in the 
space-time lattice. 

2. For each of the directions o = ±x, ±y, ±r, 
calculate the probabilities A" s . with A" a . = 
min(l,exp(-A£; s yiO), A££ = 'e% - E%. ' 

3. With uniform probability choose a direction cr. 

4. With probability A a . accept to go in the direction 
cr, and with probability 1 — A a . go to [31 
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5. Update J°, for the direction chosen and move the 
worm to the new lattice site s i+1 . 

6. If Si ^ si go to|2 

7. If Si — si go to 1 with probability po and to[3]with 
probability 1— p (p G (0, 1) and usually po = 1/2). 
We use po — 1/2 in the following. 

One advantage of this algorithm is its simplicity and the 
fact that a constructed worm is always accepted, on the 
other hand this algorithm is not directed and steps 3-4 
above are quite wasteful since in many cases the worm 
is not moved. This is avoided in the geometrical worm 
algorithm at the price of occasionally having to reject 
a complete worm. The geometrical worm algorithm, as 
described in the preceding sections, should be straight- 
forwardly applicable to the high temperature expansion 
as it was done in Ref. @ using the classical worm algo- 
rithm. We expect that this would enhance the efficiency 
of the Monte Carlo sampling. 



IV. PERFORMANCE OF THE ALGORITHMS 

Here we present results on autocorrelation times ob- 
tained with both directed and undirected algorithms. For 
the sake of brevity, we restrict ourselves to the case p = 0, 
where the critical point is known with high precision, 
and where results on autocorrelations for the undirected 
worm algorithm have already been published |6J . All the 
results presented in this section correspond to runs on 
cubic lattices of 10 7 -10 8 Monte Carlo worms for a value 
of K = 0.333, extremely close to the critical point (esti- 
mated as K c = 0.33305(5) in Ref. [(j). In principle, simu- 
lations should be performed on lattices of size Lx Lx L T , 
but here since the dynamical exponent z = 1 at fJ, = 0, 
we can set L T = L. We focus here on calculations of the 
energy E = (H) and the stiffness p defined as 



(15) 



where L is the linear size of the lattice. 

For the simulations with directed worms, we restrict 
ourselves to |J| < 3 for the tabulation of probabilities. 
Probabilities involving higher values of |J| were calcu- 
lated during the construction of the worm. Such config- 
urations were found to be exceedingly rare. 

For the case at hand, only 1% of the "scattering" ma- 
trices P Si contained diagonal elements corresponding to 
a non-zero back-tracking probability. Moreover, these 
back-tracking (bounce) processes were found to occur for 
very unlikely configurations. The acceptance rate, 1 — P e , 
is very high for both algorithms at K c (around 98% for 
undirected worms and 97% for directed worms for all 
lattice sizes). For the classical worms, all worms are ac- 
cepted due to the nature of the algorithm. However, we 
found that many proposed attempts at changing one link 
were refused (more than 60% in our simulations). 
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FIG. 2: Auto-correlation function of stiffness versus Monte 
Carlo time (defined by one worm construction - see text -) for 
L = 56 at K = 0.333 for directed (circle), undirected (dotted 
line) and classical (solid line) algorithms. 



In figure [21 we present the autocorrelation function of 
stiffness for a lattice size L = 56, for directed, undirected 
and classical worms. The autocorrelation function Co(t) 
of an observable O is defined in the standard way: 



Co{t) 



(O(t)O(O)) - {Of 
(O 2 ) - {O} 2 



(16) 



where (. . .} denotes statistical average and t is the MC 
time, measured in the number of constructed worms (ac- 
cepted or not). We clearly see in figure that the di- 
rected worm algorithm is more efficient at decorrelating 
the data than undirected and classical worms, the latter 
having the longest autocorrelation times. 

Now we define the autocorrelation time tq of an ob- 
servable O. In Ref. 0, to was defined as the greater 
time of a double-exponential fit of the autocorrelation 
function. Here we use a much simpler definition, inde- 
pendent of any fitting procedure: to is defined as the 
time where the normalized autocorrelation function de- 
crease below a threshold to ■ We can use different thresh- 
olds for different observables O. Since for small lattices 
and especially for directed worms, autocorrelation times 
are small, and since Coif) is known only for discrete val- 
ues of t, to is determined by a simple linear interpolation 
between the two times surrounding the threshold. It is 
important to note that the values of the autocorrelation 
time depends on the threshold to, but the dependence 
on lattice size of these autocorrelation times should not 
change as long as to is small enough. Error bars on to 
have been estimated by slightly changing the threshold, 
by an amount in between 2% and 5% in this work. 

Using the above mentioned determination of autocor- 
relation times, we extract autocorrelation times of the 
stiffness p and the energy E for directed, undirected and 
classical worms. The threshold was set the same for 
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FIG. 3: Auto-correlation times of the stiffness, p, for di- 
rected, undirected and classical algorithms versus lattice size 
L. Shown are the raw auto-correlations times, before rescaling 
to take into account the computational effort expended. 
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FIG. 4: Auto-correlation times of energy E for the three 
algorithms versus lattice size L. Shown are the raw auto- 
correlations times, before rescaling to take into account the 
corresponding computational effort expended. 



all algorithms when comparing the same quantity: we 
used through this work t p — 0.02 for the stiffness and 
tg = 0.05 for the energy. Scaling of these times with the 
lattice size is shown in figure|3]for stiffness and in figure^] 
for energy. It can be seen that whereas autocorrelation 
times grow approximatively linearly with lattice size for 
all algorithms, the slope is significantly smaller for the 
directed worm algorithm. 

It is clear from these results that the directed algorithm 
significantly reduces the autocorrelation times. However, 
the average size of the directed worms could be larger, 
and hence on average consume more computational time. 
For all algorithms the computational effort is linearly pro- 



FIG. 5: Mean size (to) divided by 3L 3 versus lattice size L for 
directed, undirected and classical worms. 



portional to the length of the worm. To make an honest 
comparison, we therefore have to multiply the autocor- 
relation times by the number of attempted changes per 
link, which we define as (w) / (3L 3 ) , where (w) is the mean 
worm size (the mean number of links the worm has at- 
tempted to visit), L the lattice size preceded by an irrele- 
vant factor indicating that there are 3 links per site. For 
the classical worms, the mean worm size (w) is defined 
as the total number of proposed attempts (step 4 in the 
pseudo-code presentation in section IIII C(l . In order to 
make an un-biased comparison of the three algorithms it 
is here necessary to include the updates refused during 
the construction of the classical worms in the definition 
of (w). 

As mentioned, the computational effort (the CPU 
time) is linear in (w) for all algorithms. An equivalent 
rescaling was used in Ref. 6] in order to make a fair 
comparison with the local algorithm. In Fig. is shown 
the mean worm size (w) (divided by 3L 3 ) for the three 
algorithms versus lattice size, corresponding to the aver- 
age fraction of the total number of links occupied by the 
worm. In both cases, we see that this fraction decreases 
with L. We also note that the classical worms are longer 
than in the other proposed algorithms, which will result 
in larger autocorrelation times. Directed and undirected 
worms are almost of the same size, with very slightly 
larger directed worms. The corresponding effect on the 
value of rescalcd effort (presented in the next paragraph) 
will be small when comparing autocorrelation times for 
those algorithms, however we wish to keep it present for 
a more fair analysis. 

Having discussed the behavior of (w) we now take 
into account the computational effort used to construct 
the worm by rescaling the auto-correlation times by 
(w)/(3L 3 ). We show in figure El the rescaled auto- 
correlation times for the three algorithms. We find that 
the auto-correlation times per link stay reasonably small 
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FIG. 6: Auto-correlation times of stiffness p and energy E for 
the three presented algorithms versus lattice size L. These 
auto-correlations times are rescaled auto-correlation times 
where the computational effort is taken into account. 
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FIG. 7: The probability density, P(w) for generating a worm 
occupying a fraction of w/3L 3 of the lattice, as a function of 
w/3L 3 for the directed and un-directed algorithms. Shown 
are results for a lattice of linear size L = 56 at K — K c . The 
solid lines indicate power-law fits to the data. 



for all algorithms, but the directed algorithm clearly gives 
better results, with auto-correlation times smaller by a 
factor around 4 (1.5 — 1.7) for the stiffness (energy) with 
respect to the undirected algorithm, and a factor around 
10 (4) with respect to the classical worm for the largest 
sizes. The fact that both algorithms are more efficient at 
decorrelating the stiffness than the energy seems to indi- 
cate that the worms couple more effectively to "winding 
modes" , from which the stiffness is uniquely determined, 
than to simple local modes which determine the energy. 
With the same argument, we can see that directed worms 
are more efficient at updating winding modes than undi- 
rected or classical worms. 

The actual distribution of the size of the worms gen- 
erated, P(w) is also of interest. In Fig. we show re- 
sults for the probability density, P(w) for generating a 
worm occupying a fraction of w/3L 3 of the lattice, as a 
function of w/3L 3 for the directed and un-directed algo- 
rithms. The classical worm algorithm has a distribution 
identical to the one shown for the un-directed algorithm. 
Clearly, the directed worms have a somewhat broader 
distribution but for both algorithms the distribution fol- 
lows a power-law form P(w) ~ w~ a with a ~ 1.37. The 
power-law behavior is to be expected since the simula- 
tions were performed at the critical point. Away from 
the critical point we have verified that the initial power- 
law form crosses over to an exponential behavior at large 
arguments. 

To summarize, we find that in all cases, rescaled auto- 
correlation times stay almost constant with the lattice 
size, but could also be fitted with a very small power- 
law or logarithm, showing an almost complete elimina- 
tion of critical slowing down. All in all, the main result 
of this section is that directed worms produce less cor- 
related data (smaller autocorrelation times), even if the 



scaling is good for all the three (directed, undirected and 
classical) algorithms. We also note that the geometrical 
worm algorithms perform better than the classical worm 
algorithm. 

The fact that both directed and undirected algorithms 
have almost the same scaling of the rescaled autocorre- 
lation time with L, seems also to be observed for the di- 
rected loop Quantum Monte Carlo cluster algorithms . 



V. THE CORRELATION FUNCTIONS 

A. Measurements of correlation functions with 
worm algorithms 

For the quantum rotor model, the correlation functions 
of interest have the following form : 

C(r,r / ,T,T / ) = (e i(9 "- (T )-^' (T ' )) ) , (17) 

where the O's are operators for the phase of the bosons, 
and e l6r ^ = e TH e lSr e~ TH . Due to transnational invari- 
ance, C(r,r',r, r') = C(r — r', r — r'). Physically this 
corresponds to creating a particle at (r, r) and destroying 
it at (r',r'). When this correlation function is mapped 
onto the link-current representation the creation and de- 
struction of the particle is interpreted as a particle cur- 
rent going from (r, r) to (r',r'). As is evident from the 
definition of the correlation function in Eq. <|17[) the value 
of the correlation function can not depend on the specific 
path taken from (r, r) to (r',r') as long as we take into 
account the fact that going in the x,y,r increases the 
local current, whereas going in the —x, —y, — r direction 
decreases the local current. In the link current represen- 
tation this correlation function can be written |ll| in the 



following way: 

C(r,r) = { J] 

ex P (sign(oi) (jf rtiTi) - ^.irMr,) + I ), 

(18) 

where "path" is any path on the space-lattice connect- 
ing two points a distance (r, r) apart and o~i is the di- 
rection needed to go from (ri,Tj) to (ri+i, t^+i), <r, = 
±x, ±y, ±r. When going in the direction cr, = x,y,T we 
propagate a particle and the correlation function corre- 
sponds to incrementing the corresponding link-variable 
by one. When going in the direction Oi = —x,—y,~T 
we propagate a hole in the x,y,z direction and the cor- 
relation function corresponds to decrementing the cor- 
responding link-variable by one. This is indicated in 
the expression Eq. (|18|l by sign(cTi). Furthermore, we 
only get a contribution from /i ri whenever we go in the 
r— direction and we take this into account by 5 ait ± T . If 
we define J7* y T j = ~Jf x ~iyT) with analogous defini- 
tions for the other directions we see that by incrementing 
and decrementing the link-current variables in the above 
manor ^2 a J7 TT = at all the sites between (ri,Tj) and 
( r i+ii T i+i)- The current is divergenceless at all the inter- 
mediary sites. The sites (r^r^) and (ri+i, Tj+i) will have 
non-zero divergence with J? r — 1 corresponding to 

a site where a particle is created (or a hole destroyed) . A 
site with J{r ^ = — 1 is a site where a hole is created 
(or a particle destroyed). In Fig. [HJwe show two possi- 
ble paths V a and Vb for the evaluation of the correlation 
function C(r, r). As usual, C(r, r) = C(r + L,t + L T ) 
but C(r, r) is in general not equal to C(r, — t). 

Previous work [TJ [3 have attempted to calculate the 
correlation function by evaluating the thermal expecta- 
tion value in Eq. (|18|l along a straight path from (r, r) to 
(r',r'). Although formally correct, this method fails for 
large arguments of the correlation function due to the 
fact that for a given configuration of the link-variables 
roughly only one specific path between (r, r) and (r', r') 
will yield a contribution of order 1 . 

The geometrical worm algorithm allows for a much 
more efficient way of evaluating the correlation functions. 
In essence, before the worm returns to the starting site, 
the path of the worm corresponds precisely to the cre- 
ation of a particle at site s\ and the destruction at the 
current site Sj with a current going between the two sites. 
This is precisely the Greens function that we want to cal- 
culate. More precisely we extend Eq. I|18l) to include a 
summation over all possible paths: 

<*'.'> n 

V (r 4 ,Tj) ev 

CX P |-^= fsign(o-j) (j" ruTi) - ^i-rMr,) + §J j )• 

(19) 
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FIG. 8: Two possible paths, V a and Vb, for the evaluation 
of C(r, r). When the path is going in the x,y,r direction a 
particle is propagated in the forward direction corresponding 
to an increment in the current. When the path is going in 
—x, —y, —t direction we propagate a hole in the forward di- 
rection corresponding to a decrement in the current. The solid 
circles correspond to sites where a single particle is created or 
destroyed. 

Here V is a path for the correlation function and N-p is 
the number of paths included in the sum. Since the geo- 
metrical worm algorithm generates paths between (r, t) 
and (r n ,T„) with the correct exponential factor (except 
for a multiplicative constant) it is now easy to calculate 
the correlation functions. 

Suppose that we, by using either the directed or undi- 
rected worm algorithm, have reached the equilibrium 
configuration fj,. The probability for, during the construc- 
tion of a worm starting at site si = (ri,Ti), creating a 
current j that reaches s n = (r„, r„) ^ si is given by: 

n-i » a 

P(j w ^p')=P(s 1 )l[ 1 f- (20) 

»=i Si 

for the undirected algorithm. For the directed algorithm 
we have: 




If we call the resulting state // we can calculate the prob- 
ability for, starting from //, creating an anti-current, j, 
going from s n to s\. We find for the undirected algo- 
rithm: 

2 — 

P(j; M '^ M )=P( S „)n^, (22) 

i—n 1 
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and for the directed algorithm: 



.4 



P(j\rf -> A 1 ) = P^n)^- II Pa t (si-i\si+i). (23) 



In both cases we see that 



N± TT 



P(j;»'^n) 

expj--^ ^sign(CTi) (j ( 



(r<,Ti) 



(r*,Ti) 



-5, 



(24) 



Hence, we see that for both algorithms the intermedi- 
ate states generated during the construction of the worm 
follows precisely the distribution needed apart from the 
factor N Sn / N Sl . It follows that whenever a worm reaches 
a point a distance (r, r) away from the initial point it con- 
tributes a factor of N Sl /N Sn to the correlation function of 
argument (r, r). Note that it follows from the above proof 
that all worms, even the ones that are finally rejected, 
have to be included in the calculation of the Greens func- 
tions. Per definition C*(0,0,0) = C(L,L,L T ) = 1. 
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FIG. 9: The Greens function C(x) for a system of size L s , L = 
64 at K — K c = 0.33305, /i = 0, as a function of x. The solid 
line indicates a power- law fit of the form 0.296 (a; -1 035 + (64- 

xr 1 - 035 ). 



It would be of much interest to calculate C(r, r) for 
fi 7^ using this method. Such calculations are currently 
in progress |23 ■ 



B. Results 

The above procedure is straight forward to implement. 
Suppose we want to calculate the Greens functions for a 
d+1 dimensional system with d = 2. Since the two space 
directions are equivalent by symmetry it is only necessary 
to calculate C(x,t). This is easily done by keeping track 
of the position of the worm during construction. If the 
relative position of the worm with respect to its starting 
point Si, is denoted by (x r ,y r ,T r ), when the worm has 
reached site s n , we add N Sl /N Sn to C(x r ,y r ,T r ). This 
can be done with very little computational effort and 
since an enormous amount of worms are generated during 
the simulation extremely good statistics can be obtained 
for C{x r , y r ,T r ) by averaging over the worms (which can- 
not be achieved with the local algorithm) . As mentioned, 
in order not to bias the calculation, even worms that are 
eventually rejected should be included for a correct cal- 
culation of the Greens functions. In Figure we show 
results for the Greens function as a function of x for a 
system of size L 3 ,L = 64. For this simulation the di- 
rected algorithm was used with a total number of worms 
equal to 1.5 x 10 8 . It is easy to obtain extremely small 
error bars on the Greens functions even for very large 
system sizes. For the results shown in Fig [5] /i = and 
by symmetry C(r) is identical to C{x). From scaling 
relations |2l| C(r) is expected to decay as r -( d - 2 + z +'n) 
where z is the dynamical critical exponent. With /i = 0, 
z = 1 we find C(r) ~ r _ ' 1+r " . Fitting to this form we 
find rj — 0.035(5). The obtained critical exponents are 
in excellent agreement with previous work |llj | and more 
recent high-precision estimates for the critical exponents 
of the 3d XY model HJ. 



VI. SUMMARY AND DISCUSSION 

We have proposed a directed worm algorithm for the 
quantum rotor model. This algorithm is an improve- 
ment of the "undirected" algorithm presented in It 
has been shown that by adjusting the degrees of freedom 
left in the detailed balance condition, one can construct a 
more efficient algorithm by minimizing the back-tracking 
(bounce) probability for the worm to erase itself. The 
minimal probabilities can be found by solving a linear 
programming problem subject to a few well-defined con- 
straints. A proof of detailed balance for the directed case 
has also been presented. The directed and un-directed al- 
gorithms are identical except for the fact that appropri- 
ately defined local probabilities p a for moving the worm 
through the lattice are chosen in an optimal manner 
for the directed algorithm. Hence, only a very limited 
amount of additional programming has to be done to im- 
plement the directed algorithm. 

These central ideas for this directed algorithm can be 
straightforwardly applied to directed QMC loop algo- 
rithms and one can avoid an analytical calculation 
for each new model where one wants to implement a di- 
rected algorithm. More generally speaking, we believe 
that the framework presented here could be useful for 
constructing new algorithms for other models, for exam- 
ple classical spin models |23J. 

We have shown the superiority of the directed algo- 
rithm as compared to the undirected one and to the ap- 
proach (" classical worms" ) proposed in 8] by calculating 
autocorrelation times of different observables near a crit- 
ical point. Whereas the computational gain is not as 
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drastic as when pass ing fro m a local update algorithm 
to a worm algorithm [g, [l7| , we showed that one gains a 
factor ranging from 1.5 to 10 (depending on the quantity 
and on the comparison) for the simulations considered 
here. We did not try to estimate autocorrelation expo- 
nent z for the algorithms, because in all cases, it is small 
(as can be seen in figureEJ) and it would be hard to deter- 
mine with high precision. Looking at the data, it is likely 
that values of z for all algorithms are the same or quite 
close. A logarithmic dependence of r on L, indicating 
z = 0, cannot also be excluded. 

In this paper, we have also derived an efficient way 
of measuring correlation functions during the worm con- 
structions. This feature is similar to other worm algo- 
rithms 0, 0] , but here we show, including analytical ar- 
guments, that it also works for directed worms. The situ- 
ation for directed QMC loop algorithms is less certain, 



even if some results were recently presented in Ref . |24j . 

The directed worm algorithm could be specially useful 
to study the transition for a non-commensurate value of 
the chemical potential in the pure quantum rotor model 
or for the disordered case, where very strong finite size 
effects have been identified [T3, EH |22| • 
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